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1.  Introduction 


This  report  addresses  issues  associated  with  the  modeling  and  simulation  of 
atmospheric  transport*  and  diffusion  processes.  Specifically,  currently  used 
methodologies  are  compared  and  examined  for  range  of  applicability  and  for 
appropriateness.  The  topic  should  be  of  interest  to  Department  of  Defense 
(DoD)  users  and  program  sponsors  that  rely  on  such  models  for  decision.  The 
topic  also  fulfills  a  requirement  of  the  Ballistic  Missile  Defense  Organization. 
To  be  certain  that  all  readers  have  the  same  basic  understanding,  a  substantial 
portion  of  this  report  is  concerned  with  providing  a  necessary  but  not  exhaustive 
overview  and  history  of  the  concepts  that  provide  the  foundation  for  existing 
models.  While  some  readers  may  be  uncomfortable  with  the  mathematics  and 
others  concerned  with  the  lack  of  detail,  every  effort  has  been  made  to  assure 
sufficient  coverage  that  technically  advanced  readers  are  able  to  draw  their  own 
conclusions,  while  others  achieve  some  understanding  without  being 
overwhelmed  with  mathematical  intricacies. 

To  put  things  in  the  proper  perspective,  it  is  well  to  remember:  before 
computers,  solutions  to  complex  physical  problems  could  only  reasonably  be 
treated  in  an  idealistic  parametric  manner  and  were  fraught  with  many 
simplifying  assumptions.  Still,  beyond  the  simplest  of  cases  (ones  for  which 
exact  mathematical  results  exist)  the  accuracy  of  the  result  is  obviously 
dependent  on  how  well  the  parameterization  replicates  the  physical  process. 
The  prediction  of  hazards  caused  by  the  dissemination  of  any  material  into  the 
atmosphere  is  an  extremely  complex  problem  driven  primarily  by  meteorology 
and  affected  by  source  structure,  chemical  reactions,  decay  factors,  and  human 
response.  While  each  factor  is  in  itself  a  complex  problem,  and  the  final 
solution  is  only  as  accurate  as  the  least  accurate  treatment  in  the  chain,  this 
paper  addresses  only  the  issues  associated  with  atmospheric  phenomena  that 
control  the  dispersion  of  materials,  commonly  referred  to  as  transport  and 
diffusion  (T&D). 


*Terms  in  this  report  are  referred  to  in  the  index  beginning  on  page  59. 
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Local  meteorological  conditions  are  the  primary  determinant  of  the  outcome  of 
the  T&D  processes,  and  a  complete  understanding  of  atmospheric  behavior  is 
essential  to  the  description  of  its  influence.  To  ascertain  the  effect  of 
atmosphere  on  the  mass  dissemination,  we  must  first  consider  that  forces  of 
much  different  temporal  and  spatial  scales  drive  these  two  separate  yet  coupled 
physical  processes.  The  scale  and  coupling  are,  in  turn,  dependent  on  space  and 
time  considerations.  [1]  Using  the  standard  notion  of  transport  (the  motion  of 
the  center  of  mass  of  some  defined  entity)  does  not  relieve  the  necessity  to 
examine  the  temporal  and  spatial  variability  that  leads  to  a  three  dimensional 
smearing  about  that  center.  Traditionally,  we  treat  all  motions  not  associated 
with  the  motion  of  the  centroid  as  diffusion,  and  we  further  differentiate 
molecular  or  microdiffusion  from  macroscopic  turbulent  diffusion.  This 
separation  is  natural  because  molecular  diffusion  is  a  property  of  the 
transporting  media  while  turbulence  is  a  property  of  the  flow.  The  effects  of 
turbulent  diffusion  are  generally  orders  of  magnitude  larger  than  the  effects  of 
molecular  diffusion,  a  fact  many  researchers  use  to  justify  ignoring  the  later. 
As  will  be  seen,  the  common  approach  is  to  mathematically  treat  turbulent 
diffusion  in  a  manner  similar  to  molecular  diffusion.  This  has  led  to  many 
physical  constructs  that  attempt  to  put  order  into  a  process  that  is  inherently 
unordered.  A  pervasive  notion  is  that  turbulence  is  a  random  process;  therefore, 
it  eludes  a  rigorous  mathematical  description  necessitating  an  empirical  or 
statistical  treatment.  To  gain  an  appreciation  for  the  various  methods  for  the 
treatment  of  this  phenomenon,  it  may  be  best  to  review  the  mechanisms  that 
introduce  turbulence  into  the  atmosphere  and  the  various  methods  used  to 
characterized  it  and  its  effect. 


2.  Atmospheric  Motions 


In  the  absence  of  differential  energy  effects  (surface  heating  and  density 
variations),  a  very  simple  picture  of  a  fluid  surrounding  a  spinning  globe  can 
be  described  analogous  to  well-studied  viscous  flow  over  a  flat  plate.  For  a 
smooth,  uniform  surface  and  relatively  slow,  undisturbed  flow,  two  regions  can 
be  defined.  The  lower  level  is  the  region  directly  influenced  by  the  frictional 
drag  (no-slip  condition)  at  the  fluid-solid  interface  and  is  characterized  by  an 
increasing  velocity  for  some  distance  from  the  surface.  The  flow  of  the  upper 
region  is  independent  of  frictional  effects  resulting  in  a  constant  velocity  with 
distance  from  the  surface.  In  both  regions,  the  flow  is  typified  as  laminar, 
which  can  be  visualized  by  the  introduction  of  a  steady  source  of  dye.  In 
laminar  flow,  the  dye  traces  a  thin  ribbon  or  streamline  in  the  direction  of  the 
flow,  horizontal  and  vertical  spreading  of  the  dye  is  caused  by  molecular 
diffusion  and  the  transfer  of  momentum  in  the  region  influenced  by  friction. 

Irregular  surfaces  (hills,  grasslands,  and  wooded  areas)  disturb  the  purely 
laminar  aspects  of  the  flow  in  two  ways.  Experiments  by  Osborne  Reynolds 
in  the  late  nineteenth  century  demonstrated  that  flow  over  a  rough  surface  could 
be  characterized  by  a  nondimensional  number  used  to  indicate  various 
transitions  in  flow  characteristics.  The  Reynolds  number  (Re)  results  from  a 
ratio  between  the  elements  that  contribute  to  and  suppress  a  disturbance  in  the 
flow.  Roughness  and  speed  lead  to  an  unsteadiness  while  viscosity  tends  to 
dampen  it,  leading  to  Re  =  lU/v,  where  1  characterizes  the  length  (height)  of  the 
disturbing  element(s),  U  is  the  mean  flow  speed,  and  v  the  kinematic  viscosity 
of  the  fluid.  Reynolds  experiments  showed  that  for  very  small  Re,  the  flow 
remained  laminar  in  the  sense  that  streamlines  remained  parallel.  Returning  to 
the  flow  over  a  plate,  in  this  case  with  surface  protuberances,  the  maximum 
height  at  which  flow  velocity  continues  to  change  is  simply  higher.  Two 
distinct  flow  regimes  are  still  evident. 

The  experiments  further  showed  that  as  Re  increased,  circulating  eddys  formed 
on  the  downstream  side  of  an  obstruction,  and  at  even  larger  Re,  the  eddys  shed 
into  the  laminar  flow  disturbing  its  uniform  characteristics  downstream  of  the 
obstacle.  Eddy  size  and  the  periodicity  of  shedding  are  directly  correlated 
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to  Re.  Extending  these  results  to  atmospheric  flow  in  proximity  to  the  very 
nonuniform  surface  of  the  earth  and  in  the  absence  of  any  thermal  effects,  the 
production  of  a  very  large  range  of  eddy  sizes  is  possible.  This  action  is  one 
of  the  primary  sources  of  mechanical  turbulence  in  the  atmosphere,  providing 
the  rationale  for  characterizing  turbulence  as  eddy  turbulence  and  resulting 
dispersion  as  eddy  diffusivity.  As  is  discussed,  turbulence  is  not  as  well 
organized  as  this  view  may  suggest;  it  does,  however,  provide  a  basis  for  the 
definition  of  turbulence  regimes  in  the  atmosphere.  The  complete  picture 
requires  that  we  also  consider  thermal  eddy  sources,  and  the  effects  by  large  air 
masses  distributed  throughout  our  atmosphere. 

The  simplified  notion  described  is  complicated  first  by  temperature-density 
variations  and  by  differential  heating  of  the  surface  of  the  earth. 
Temperature-density  variations  are  manifested  in  a  tendency  for  the  atmosphere 
to  organize  itself  into  large,  identifiable  air  masses.  The  boundaries  between 
air  masses  or  weather  systems  can  result  in  velocity  shears  and  strong  gradients 
providing  additional  sources  of  turbulent  energy.  Each  air  mass  type,  exclusive 
of  its  boundaries,  has  its  own  turbulence  structure.  Superimposed  on  the  large 
scale  features  are  the  effects  of  the  local  surface  topography  and  heating. 
Surface  heating  creates  small  convective  or  thermal  eddys  that  provide  a 
mechanism  for  the  vertical  transport  of  turbulence.  In  the  absence  of  the 
thermal  effects,  turbulent  eddys  caused  by  surface  irregularities  would  remain 
confined  to  a  very  thin  layer  close  to  the  surface.  In  the  real  atmosphere, 
vertical  convection  caused  by  surface  heating  transports  the  eddys  to  higher 
levels. 

Based  on  these  considerations  a  reasonable  separation  of  the  atmosphere  into 
regimes  having  very  different  turbulence  structures  is  possible.  Meteorologists 
identify  a  surface  boundary  layer  (SBL)  in  which  surface  features,  friction, 
viscosity,  and  temperature  gradients*  are  the  dominant  influence  on  the  motion. 
Above  this  layer  is  a  transition  zone  or  planetary  boundary  layer  (PBL)  in 
which  density  gradients,  residual  friction  (caused  by  convective  transport)  are 


‘Strongly  heated  surfaces  create  temperature  gradients  near  the  surface  much  larger  than 
those  found  elsewhere  in  the  atmosphere.  In  some  cases,  the  gradients  lead  to  isolated 
rivers  of  vertical  flow  that  cause  mirages. 
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the  dominant  flow  determinants.  Finally,  far  removed  from  the  surface  we  have 
a  region  that  can  be  treated  as  inviscid  because  the  governing  forces  are 
primarily  pressure  gradients  and  the  coriolis  effect  (caused  by  rotation  of  the 
earth).  The  highest  layer  is  generally  termed  the  free  atmosphere,  free  in  the 
sense  that  motions  are  all  but  independent  of  surface  characteristics  and 
influences. 


The  thickness  of  each  layer  is  highly  dependent  on  prevailing  conditions  (terrain 
roughness,  wind  speed,  temperature,  time  of  day/year,  cloud  cover,  etc.),  but  in 
general,  the  SBL  extends  50  to  100  m  above  the  surface  and  the  planetary  layer 
extends  500  to  1000  m.  Using  standard  fluid  dynamic  analyses,  [2]  we  can 
express  these  definitions  in  a  more  rigorous  manner.  Accelerations  at  any  point 
in  an  incompressible  fluid  are  described,  in  Cartesian  coordinates,  by 
Newton’s  second  law. 


./du  du  du  du x  ,&u 

p(—  +  u —  +  v —  +  w — )  =  |i( -  + 

dt  dx  dy  dz  dx2  dy2 


cPu  dP  _ 

—  +  TT)  '  T"  +  PS 

dz 2  dx 


(1) 


where 


u,  v,  and  w 

P 

P 


=  the  flow  velocities  in  the  respective  x,  y,  and  z  Cartesian 
coordinates 
=  the  fluid  density 
=  the  dynamic  viscosity 
=  the  pressure 

=  other  external  forces  per  unit  mass. 


Similar  equations  for  the  y  and  z  components  yield  three  equations  that  form 
the  Navier-Stokes  equations  of  motion.  The  term  on  the  left-hand  side  is  the 
inertial  or  acceleration  term  (the  ma  of  F  =  ma).  On  the  right  are,  the  viscous, 
pressure,  and  additional  (coriolis  force  and  gravity)  force  terms,  respectively. 
It  is  the  viscous  force,  together  with  the  physical  reality  that  at  a  surface,  the 
relative  velocity  must  vanish  (no-slip  condition),  which  accounts  for  resulting 


‘For  Reynolds  numbers  associated  with  atmospheric  motions,  the  simplifying 
approximation  of  an  incompressible  atmosphere  has  been  shown  to  be  valid. 
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frictional  and  rotational  forces  that  lead  to  eddys  and  turbulence. 
Mathematically,  the  viscous  term  can  be  separated  into  two  terms  reflecting 
frictional  and  rotational  forces.  Therefore,  the  amount  of  energy  existing  in  the 
turbulent  eddys,  created  by  the  rotational  forces,  can  be  tracked.  The  inertial 
term  is  nonlinear,  which  makes  universal  solutions  extremely  difficult  except 
in  the  simplest  cases.  [3]  Separating  the  atmosphere  into  regimes  based  on 
importance  of  the  various  force  terms  permits  convenient  generalizations  for  a 
layer.  Unfortunately,  success  in  applying  generalizations  in  one  domain  does 
not  guarantee  the  similar  applicability  in  other  regimes. 


In  the  SBL,  the  temperature  gradient  and  viscous  terms  in  equation  (1)  can  be 
shown  to  be  much  larger  than  either  the  pressure  or  coriolis*  terms.  The 
manner  in  which  the  temperature  gradient  is  included  will  be  shown  in  the 
section  on  diffusion.  As  distance  from  the  surface  increases,  the  viscous  and 
temperature  gradient  terms  decrease  to  a  point  at  which  the  order  of  all  terms 
is  the  same.  The  process  can  be  visualized  as  follows.  Small  volumes  of  air 
immediately  adjacent  to  the  surface  are  heated  by  conduction,  and  rise  similarly 
to  bubbles  in  boiling  liquid.  Remaining  intact  and  rising  through  ever 
decreasing  external  pressure  forces,  parcels  cool  through  expansion;  and  moving 
slower  than  the  surrounding  air  act  as  obstacles  in  the  flow,  effectively 
extending  the  influence  of  the  surface  through  a  larger  than  otherwise  layer. 
The  rising  volumes  constitute  an  eddy  and  can  encompass  eddys  formed  by  the 
horizontal  flow  near  the  surface.  The  height  M  to  which  these  thermally  driven 
eddys  rise  (transporting  horizontal  momentum,  heat,  and  turbulence)  depends 
on  the  initial  impulse  and  the  diminishing  buoyancy  force  resulting  from  density 
differences.  M,  the  effective  top  of  the  PBL,  referred  to  as  height  of  the  mixed 
layer,  contains  the  planetary  and  surface  layers.  The  ability  of  the  atmosphere 
to  produce  such  a  mixed  layer  can  be  described  through  the  ratio  between 
competing  convective  and  advective  forces.  The  Richardson  s  number  R;, 


*The  coriolis  force  is  of  the  order  cov,  where  co  is  the  angular  velocity  of  the  earth  s 
rotation,  and  v  is  the  flow  velocity. 


8 


where 


g  =  gravitational  acceleration 

r  =  the  appropriate  adiabatic  temperature  lapse  rate  (see  footnote  on 
page  20). 


Comparing  actual  temperature  lapse  rates  to  the  adiabatic  rate  is  a  standard 
meteorological  means  of  determining  the  vertical  stability  of  an  atmospheric 
layer.  Examination  of  equation  (2)  indicates  that  single  point  measurements  of 
T  and  u  are  inadequate  descriptors  of  a  mixed  layer.  To  correctly  characterize 
the  mixing  potential,  realistic  gradients  of  these  quantities  are  required.  For 
Ri  >  1,  turbulence  is  suppressed  and  a  tendency  to  laminar  flow  is  achieved;  for 
Ri  <  1,  the  flow  tends  to  remain  turbulent  and  the  mixed  layer  tends  to  grow 
in  thickness. 


3.  Atmospheric  Turbulence 


Section  2  suggests  that  mechanisms  introduce  turbulence  into  the  atmosphere, 
and  indicates  that  the  random  nature  of  turbulence  lends  it  intractable.  If  the 
processes  are  completely  random,  a  less  rigorous  but  pliant  description  can  be 
applied  using  statistical  methods.  Before  attempting  this,  it  is  good  to  isolate 
the  tractable.  Following  Reynolds,  this  can  be  accomplished  through  the 
definition  of  any  variable  term  through  the  superposition  of  its  mean  (a)  and 
fluctuating  (a')  components,  where: 

t  +  .5  T 

J  a(t)dt 

a'  =  a  -  a;  a  =  -  - - .  (3) 

t  *  .ST 

f  dt 

t  -  ,5T 


Applying  this  notion  to  equation  (1),  it  is  possible  to  separate  the  flow  into  its 
mean  and  fluctuating  components.  Clearly,  the  viability  of  this  approach  is 
dependent  on  choosing  a  period  T  such  that  a  represents  the  mean 
characteristics  of  the  general  flow.  T  must  be  sufficiently  long  to  include  an 
adequate  number  of  low-frequency  fluctuations,  while  an  overly  long  period 
may  effectively  conceal  important  features  of  the  higher  frequency  variability. 


Imposing  a  requirement  of  steady  mean  flow,  defined  by 


?-! 


t  +  .5  T 


J  (u  -  u)dt  =  u  -  u  =  0 


t  -  .5T 


leads  to  u'  =  v'  =  w'  =  0.  However,  fluctuation  correlations  (eddy  velocities 
represented  by  u'2,  u'w',  etc.)  do  not  necessarily  vanish  and  lead  to  a  problem 
in  the  solution  of  the  Navier-Stokes  equation  known  as  closure.  The  closure 
problem  is  a  mathematical  result  of  having  insufficient  information  (equations) 
to  match  the  number  of  unknown  quantities.  Thus,  even  under  ideal  conditions 
of  steady  flow,  energy  can  exist  in  the  turbulent  structure  and  contribute  to 
diffusion.  If  the  mean  flow  is  not  steady,  other  restrictive  conditions  must  be 


11 


applied.  In  all  cases,  equation  (3)  must  be  satisfied  exactly  or  to  a  high  degree 
of  approximation  before  this  construct  is  useful.  As  appealing  as  this  approach 
may  first  appear  (separation  into  a  well  behaved  and  solvable  equation  and  a 
coupled  but  unsolvable  transient  equation),  it  still  does  not  permit  a  rigorous 
definition  of  the  effects  of  turbulent  diffusion  for  one  must  have  a  priori 
knowledge  of  or  a  valid  method  for  estimating  the  transient  eddy  velocities  or 
at  least  the  amount  of  energy  contained  in  the  turbulence.  Given  the 
complexities  of  atmospheric  motions  compounded  by  weather,  terrain,  and 
diurnal  and  seasonal  variability,  the  former  is  certainly  formidable.  A  number 
of  methods,  including  statistical  approximations,  for  estimating  eddy  velocities 
or  their  effects  have  been  suggested.  These  will  be  reviewed;  but  first,  it  must 
be  noted  that  the  analyses  of  atmospheric  turbulence  are  complicated  by  the 
range  of  temporal  and  spatial  scales  over  which  the  phenomenon  occurs  in  the 
atmosphere.  Further  it  must  be  realized  that  the  regimes,  as  defined,  have 
completely  different  spectral  forms.  As  noted  by  Pasquill,  [4]  the  statistical 
properties  depend  considerably  on  the  sampling  duration,  and  often, 
autocorrelations  (e.g.,  u'2)  and  variances  do  not  satisfy  boundary  conditions. 
This  is  caused  by  variability  in  the  sources  of  turbulence  such  that  the 
atmosphere,  particularly  in  the  boundary  layer,  is  not  steady  or  stationary.  Both 
conditions  are  required  for  strict  application  of  statistical  theory,  and  more 
strictly  for  the  generalization  of  analytic  descriptions.  A  further  restriction, 
homogeneity,  is  in  force  when  the  theory  is  applied  to  the  diffusion  of  volume 
sources.  Despite  these  cautions,  statistical  approaches  are  appealing  in 
theoretical  applications  when  energy  conservation  can  be  applied  or  its 
principals  are  used  to  obtain  statistical  quantities  (e.g.,  u'2,  u'w').  One 
technique  has  been  successfully  applied  in  recent  diffusion  modeling  using  the 
Navier-Stokes  equations. 

The  atmosphere  can  be  considered  an  enormous  energy  conversion  machine, 
and  it  was  G.  I.  Taylor  [4]  who  first  suggested  that  eddys  in  the  atmosphere 
provide  the  mechanism  for  the  transfer  of  the  mechanical  energy,  contained  in 
the  large  scale  flow,  to  thermal  energy.  He  postulated  that  the  process  consisted 
of  a  cascade  of  energy  from  the  large  eddys  to  sufficiently  small  motions  that 
are  indistinguishable  from  molecular  motion  induced  by  thermal  properties. 
Kolmogorov’s  [6]  theoretical  extension  of  Taylor’s  [5]  work  provided  a 
statistically  grounded  mathematical  description  of  Taylor’s  [5]  evolutionary 


process.  The  success  of  this  idea  and  subsequent  work  is  realized  in  its  ability 
to  provide  verifiable  correlations  between  the  mean  flow  and  the  turbulence 
structure.  A  physical  realization  of  this  process  can  be  gained  by  the 
observation  that  large  atmospheric  eddys  produced  by  intermittent  sources  are 
nonpersistent,  and  are  illustrated  by  considering  a  thermally  isolated  system 
containing  a  fluid.  Any  instantaneous  agitation  of  the  fluid  that  results  in  a 
large-scale  motion,  such  as  an  abrupt  movement  of  a  wall  or  perhaps  a  stirring 
action  that  creates  a  whirlpool,  will  eventually  be  followed  by  the  (slightly 
warmer)  fluid  coming,  once  again,  to  rest. 
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4.  Diffusion 


Given  the  random  nature  of  turbulent  diffusion,  the  question  is:  "Is  there  a 
method  available  that  can  mathematically  replicate  the  physical  process?"  The 
hypothesis  by  Taylor  [5]  provides  one  link,  albeit  microscale  and  statistical  in 
nature,  between  the  thermal  and  mechanical  properties  of  the  atmosphere.  The 
atmosphere  is  a  gas;  therefore,  the  pressure  term  in  the  equation  of  motion 
equation  (1),  coupled  with  the  gas  law  (P  =  pRT)  affords  the  macroscopic 
relation  between  the  two  manifestations  of  energy.  Rendering  an  effective  and 
useful  characterization  of  the  atmosphere  at  either  scale  lies  in  our  ability  to 
provide  a  description  that  bridges  the  gap  between  them.  Much  of  the  early 
work  in  modeling  diffusion  processes  ignored  this  important  factor.  The 
conservative  nature  of  this  physical  process,  specifically  energy  conservation, 
together  with  the  many  observations  of  the  correlation  between  the  state  of  the 
mean  flow  and  the  state  of  the  turbulence  [1,4]  yields  an  appealing  method  for 
this  connection.  A  satisfactory  employment  of  this  methodology  requires  an 
adequate  scheme  to  maintain  a  close  inventory  on  energy  conversions  in  the 
domain  of  interest,  particularly  at  the  boundaries,  and  in  the  presence  of  cloud 
formation  and  water  vapor  transport.  In  the  absence  of  water  vapor 
considerations,  thermal  energy  changes  can  be  expressed  in  terms  of  variations 
in  temperature  T,  a  measure  of  heat  (energy),  by  summing  the  source  and  sinks 
of  heat  energy  in  a  volume: 


„  8T  si  r  dT  dT  dr,  XT  ir&T  &T  d2^ 
P  vdf  P  vL  djc  dy  dz  dx2  dy2  dz2 


(5) 


where 


p  =  the  media  density 
Cv  =  a  constant  specific  heat 
k  =  thermal  conductivity. 


The  first  term  on  the  right  is  the  transport  (advection)  term;  N  represents  terms 
for  energy  sources/sinks  (viscous  dissipation);  the  last  term  is  the  conduction 
term.  The  conduction  term  requires  special  attention,  and  to  do  this,  it  is 
instructive  to  derive  equation  (5). 
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The  flux  of  heat  per  unit  time  across  a  unit  area  Q  =  UpCvT  (U  is  the  mean 
velocity  in  the  direction  of  the  flux),  can  be  expressed  by 


=  -(k‘Vr){. 


(6) 


Equation  (6),  where  the  subscript  indicates  summation  in  the  three  coordinate 
directions,  implies  that  the  flux  of  heat  is  directly  proportional  to  the  gradient 
of  temperature,  with  k  representing  the  proportionality  factor.  If  the  conducting 
material  is  isotropic  (conductivity  independent  of  direction),  the  equation  can 
be  simplified  by  the  dropping  of  the  subscripts.  A  similar  equation  can  be 
derived  for  the  flux  of  a  nonuniformly  distributed  material.  The  flux  of 
material  per  unit  area  per  unit  time  J  is 

J,  -  -D, =  -(C-VX),  (7) 

OX; 


where 

X  =  the  concentration  of  the  material. 

Analogous  to  the  conduction  of  heat  in  a  solid,  D  represents  a  proportionality 
between  the  flux  of  a  material  by  self-diffusion  and  the  spatial  gradient  of  the 
material.  In  both  cases  the  flux  of  the  physical  quantity  is  an  attempt  by  the 
system  to  achieve  equilibrium  and  is  assumed  to  be  driven  by  a  force 
proportional  to  the  gradient  of  the  quantity  not  in  equilibrium.  Defining  the 
thermal  diffusivity  k  =  k/(pCv)  and  applying  the  total  derivative 
[D/Dt  =  5/5t  +  V  o  ]  (for  a  more  detailed  exposition  see  Sutton  [1]) 
equation  (6)  becomes 

—  +  V-(t/7)  =  -  V-(kV7).  (8) 

dt 

This  is  equivalent  to  equation  (5)  with  one  major  exception;  the  thermal 
conductivity  in  equation  (5)  is  outside  of  the  spatial  derivative  operator, 
implying  an  isotropic  assumption.  Equation  (7)  is  the  more  general  form  of  the 
equation  for  the  transport  and  diffusion  of  heat.  Likewise,  the  general  equation 
for  selfdiffusion  (molecular)  of  material  in  a  gas  is 


.  —  +  V-  (t/X)  =  -V-0VX), 
0f 


(9) 


and  moving  the  diffusivity  term  outside  of  the  gradient  operator  is  only  valid 
for  conditions  of  isotropy  (Fickian  diffusion).  This  is  a  simplification  that 
limits  many  of  the  transport  and  diffusion  models  in  current  use.  The  critical 
error  is  not  in  the  application  to  molecular  diffusion,  which  would  be  serious 
enough  if  used  indiscriminately,  but  in  the  extension  of  this  concept  to  the 
treatment  of  turbulent  diffusion. 

The  discussion  of  diffusion,  so  far,  is  strictly  applicable  only  to  molecular 
diffusion.  In  this  context,  the  concept  of  proportionality  between  the  force  and 
its  effect  is  physically  reasonable  and  well  founded.  Broadening  the  definition 
of  D  to  allow  it  to  represent  turbulent  diffusion,  which,  as  will  be  shown,  is  the 
basis  of  many  models,  jeopardizes  the  general  validity  of  equation  (9).  The 
scientists  who  first  suggested  the  broader  interpretation  were  cautious  in  its 
usage.  When  D  represents  molecular  diffusion,  a  linear  relationship  between 
X(x,y,z,t)  and  V%  is  reasonable  and  has  been  shown  to  be  valid.  However,  a 
simple  proportionality,  when  D  represents  turbulent  diffusion,  is  questionable; 
the  linkage  between  the  rate  of  change  and  the  gradient  has  vanished.  It  is  no 
longer  the  gradient  force  that  is  responsible  for  the  change,  rather  it  is  turbulent 
diffusion;  a  separate  and  independent  physical  process.  Although  this  is  not  as 
serious  as  might  first  be  expected*,  it  is  probably  the  singular  fact  that  accounts 
for  the  inability  of  early  attempts,  using  a  molecular  diffusion  analogy,  to 
model  the  effects  of  turbulent  diffusion.  These  failures  led  to  searches  for  the 
universal  factor.  The  constructs  that  resulted  from  these  attempts  are  so 
ubiquitous  that  they  are  often  treated  with  the  reverence  that  should  be  reserved 
for  real  physical  quantities. 


Tt  can  be  shown  that  equation  (9)  may  be  considered  valid.  However  D  is  no  longer  a 
simple  proportionality  constant  but  a  modifying  (and  complicated)  tensor-based  expression 
(Kjj)  of  a  physical  force. 
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5.  Modeling  Diffusion 


While  the  earliest  approaches  to  modeling  turbulent  diffusion  were  less  than 
successful,  a  number  of  useful  ideas  emerged  from  the  work  of  a  collection  of 
well-known  scientists.  The  constructs  developed  introduced  ideas  and 
terminology  (mixing  length  and  friction  velocity)  that  have  outlived  the 
theories,  necessitating  a  review  to  put  them  in  the  proper  perspective.  An  early 
approximation  of  turbulent  diffusion  borrowed  an  approach  used  in  analyses  of 
turbulent  vertical  transport  of  horizontal  momentum.  Defining  an  eddy 
diffusivity  K  analogous  to  the  definition  of  molecular  diffusivity,  and  replacing 
T  in  equation  (6)  by  r\,  where  r|  is  the  concentration  of  any  physical  species 
(momentum,  or  suspended  material);  Q  represents  the  flux  of  that  quantity. 
This  provided  a  basis  for  analyses  that  proceeded  on  strictly  heuristic 
arguments.  The  flux  Q'  =  KVr|',  caused  by  the  fluctuating  components,  was 
defined  as  Q'  =  (fj'u^fj'  in  which  the  prime  has  its  usual  meaning  and  u( 
indicates  the  velocity  in  the  Xj  direction.  Further,  (with  k  in  equation  (6) 
replaced  by  an  eddy  diffusivity  K)  it  was  argued  that  tj'  was  proportional  to 
(and  in  the  direction  of)  the  material  gradient  leading  to  K  =  u'/L  Later 
Prandtl,  [7]  using  dimensional  arguments*,  identified  A,  as  a  mixing  length 
analogous  to  the  mean  free  path  in  molecular  diffusion  theory.  [8]  While  this 
seemingly  added  physical  validity  t‘o  the  approach,  now  known  as  K-theory,  no 
tangible,  measurable  entity  corresponding  to  X  exists,  and  K  continued  to  be 
defined  on  empirical  grounds.  The  introduction  of  friction  velocity  u*  was  a 
further  attempt  at  defining  a  physically  reasonable  and  measurable  quantity  that 
could  be  used  to  characterize  turbulence. 

It  was  demonstrated  that  u„,  defined  as  the  square  root  of  the  ratio  between 
horizontal  shear  stress  and  material  density,  was  relatively  constant  with  height 
very  near  the  ground,. and  could  be  related  to  the  velocity  correlation  u'w\ 
Careful  studies  of  flows  over  flat  plates  showed  that  an  essentially  logarithmic 
relationship  characterized  the  ratio  u/u..  The  promise  offered  by  these  ideas, 
unfortunately,  is  limited  to  flat  terrain  and  a  narrow  set  of  atmospheric 


’Molecular  diffusivity  has  the  dimensions  of  l2/t,  if  K  is  to  be  treated  analogous  to 
diffusivity,  it  should  have  similar  dimensions. 
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conditions.  This  is  evidenced  by  the  numerous  attempts  and  equations  available 
to  obtain  u*.  However,  in  the  theory  of  momentum  transport,  studies  along  this 
line  continued,  and  it  was  shown  for  adiabatic*  conditions  that  the  eddy 
diffusion  coefficient  could  be  expressed  as  KM  =  u«Jkz.  It  was  then  argued  that 
such  a  formulation  was  equally  applicable  for  any  physical  entity.  The 
persistent  occurrences  of  nonadiabatic  boundary-layer  conditions  confounds  the 
application  of  K-theory  in  the  real  atmosphere,  which  was  derived  from 
experimental  observations  under  ideal  neutrally  buoyant  conditions.  For 
nonstable  conditions  a  number  of  functions  KM(u*,  RJ  have  been  proposed.  An 
alternative  approach  to  the  definition  of  X  (mixing  length)  was  introduced  by 
von  Karman  [9]  (also  see  Tatarski  [10]).  This  led  to  an  expression  for  u.  in 
terms  of  another  parameter,  roughness  length  (Zq). 

Numerous  attempts  at  resolving  the  inadequacies  of  K-theory  have  been 
attempted,  each  introducing  a  new  parameter  to  correct  the  deficiencies 
introduced  by  attempts  to  generalize.  The  fault  in  all  the  attempts  was  not  the 
ineptitude  of  the  scientists  but  the  difficulty  in  describing  the  behavior  of  a 
chaotic  atmosphere  in  terms  of  well-behaved  relationships  designed  for  simple 
cases  and  that  worked  well  for  relatively  quiescent  periods  (night,  flat  terrain, 
neutral  buoyancy).  The  nature  of  turbulence  and  the  inexact  method  required 
to  assign  values  to  z0  and  u*  imposes  considerable  problems  in  applications  of 
this  theory.  The  inexactness  results  from  the  fact  that  these  parameters  were 
constructed  to  represent  artificial  concepts  rather  than  measurable  physical 
quantities.  The  inclusion  of  the  Richardson  number,  which  acknowledges  the 
importance  of  the  vertical  temperature  gradient  in  determining  the  overall 
diffusion  effect,  must  be  noted;  for  in  a  highly  unstable  atmosphere,  vertical 
diffusion  is  likely  to  be  the  dominant  factor  in  the  dilution  of  the  material. 
Conversely,  an  overly  stable  condition  will  dampen  any  tendency  for  the 
material  to  diffuse  vertically.  However,  a  simple  appeal  to,  or  over  emphasis 
on,  atmospheric  stability  conditions  as  anything  more  than  a  rule  of  thumb 
masks  the  role  of  turbulence  as  the  primary  influence  on  material  dispersion. 


‘Adiabatic  conditions,  also  referred  to  as  neutral  buoyancy  conditions,  are  achieved  when 
the  vertical  temperature  lapse  rate  in  any  layer  is  equal  to  the  rate  of  change  of 
temperature,  Td  (=  5T/Sz),  that  a  dry  parcel  of  air  would  experience  if  it  were  to  rise 
without  exchanging  heat  with  its  surroundings.  Td  «  9.76  °C/km. 
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It  is  the  random  nature  of  turbulence,  coupled  with  a  questionable  correlation 
to  molecular  diffusion,  that  threatened  the  further  generalizations  of  the 
approaches  so  far  discussed,  and  led  researchers  to  return  to  the  consideration 
of  a  statistical  approach. 

Atmospheric  processes  are  naturally  Lagrangian*  and  only  weakly  correlated 
to  past  interactions.  In  contrast,  analytic  solutions  to  the  physical  equations  and 
standard  measurements  of  physical  quantities  are  Euleriant  (time  histories  of 
events  at  a  fixed  location).  Thus,  Eulerian  statistics  are  time  correlations,  while 
Lagrangian  correlations  are  spatial  in  nature.  Practically,  interpretive  results  are 
generally  portrayed  in  Eulerian  space,  while  formal  expressions  of  turbulence 
are  Lagrangian.  Therefore,  the  efficacy  of  necessary  time-space  transformations 
must  be  carefully  analyzed.  Pasquill  [4]  argues  that  errors  introduced  by  this 
transformation  are  small,  but  his  arguments  are  based  on  carefully  controlled 
experiments  in  homogeneous-isotropic-stationary  (h-i-s)  turbulence  and  limited 
observations.  He  also  points  to  the  difficulty  in  making  measurements  of  the 
Lagrangian  time-scale  and  the  wide  variability  in  the  limited  data  available 
(which  may  indicate  that  such  time-space  correlations  are  not  well  in  hand). 
Taylor’s  hypothesis,  [11]  a  utilitarian  approach  that  attempts  to  avoid  the 
time-space  problem,  provides  the  basis  for  many  turbulent  diffusion  models. 
The  tenet  of  the  hypothesis,  an  assumption  that  turbulence  at  any  point  is 
unaffected  by  a  translation  in  space  by  the  mean  wind,  can  be  shown  to  be 


*The  Lagrangian  correlation  relating  variations  with  respect  to  space  implies  that  the  effect 
on  a  particle  is  determined  by  its  physical  location  at  a  given  time  and  can  be  expressed 
by 


m  = 


u'(x, tx)  -  u'(x  +  dx, tj) 


where  5x  =  u'(x,  t0)  8t. 

Eulerian  correlations  are  based  on  measurements  at  the  same  location  at  different  times 

R(JC)  =  Ax^u'jx^) 
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tenuous  in  even  the  most  benign  conditions  for  which  the  hypothesis  is 
supposed  to  be  particularly  relevant*. 

The  Eulerian  approach  seeks  solutions  of  equation  (9)  with  possible  terms  for 
sources  and  sinks  included.  The  strength  of  this  approach  is  its  capability  to 
address  nonlinear  reactions  (sources  and  sinks)  and  wide  applicability  when  a 
suitable  function  for  KM  is  obtained.  The  Lagrangian  approach  relies  on  the 
ability  to  define  a  probability  density  function  (PDF)  i|/(Xj,t)  where  i \)  is  the 
probability  that  a  particle  is  in  a  given  volume,  or  that  it  is  at  a  location, 
Xj  <  Xj  <  Xj+dXj.  Then,  if  the  concentration  c(x,t)  is  composed  of  m-number  of 
particles, 

M  00  OO  00 

ttl 

<c(xvx2,xyt)>  =  Y  ff  / Q(*i, f/) yk(x!,t')dxl  0°) 

k  —  \  —00  —  00  —00 


where 

Q  =  the  transition  PDF  that  a  particle  at  x'  at  t'  will  be  at  x  at  time  t. 

o  indicates  an  ensemble  average  (a  statistically  averaged  result).  The  form  of 
and  ability  to  evaluate  the  transition  PDF  is  crucial  to  obtaining  Lagrangian 
solutions,  and  is  the  foremost  limitation  of  this  technique. 
Monin  and  Yaglom  [12]  argued  that  empirical  data  indicate  that  Q  is  a 
multidimensional  normal  distribution. 

Utilizing  this  assumption  of  normally  distributed1  turbulent  energy  and  further 
restricting  turbulence  to  be  homogeneous,  stationary,  and  isotropic  provide 
sufficient  simplification  that  analytic  solutions  to  equation  (10)  are  obtainable 
for  specific  cases.  Further  details  of  this  approach  are  in  the  section  on  model 


*This  topic  is  polemic  between  atmospheric  physicists.  However,  the  controversy  can  be 
lessened  when  viewed  in  terms  of  comparison  of  Eulerian  and  Lagrangian  time  scales. 

formal  distributions  are  sometimes  referred  to  as  Gaussian  distributions.  The 
requirement  that  turbulent  energy  be  normally  distributed  leads  to  Gaussian  puff/plume 
solutions. 
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descriptions.  The  obvious  questions  are:  How  representative  are  results  based 
on  these  restrictive  assumptions?  Can  a  PDF,  Q(X,t|X't'),  be  determined  for  the 
more  realistic  atmospheric  boundary-layer  case  of  nonstationary, 
nonhomogeneous,  and  nonisotropic  turbulence?  Extensive  experimental 
measurements  could  provide  these  PDFs,  but  the  validity  (adjustment  to  the 
shape  factors  for  the  Gaussian  distributions)  are  limited  to  those  specific 
situations.  Such  applications  are,  at  best,  diagnostic,  and  attempts  to  reproduce 
anything  other  than  trends  are  irresolute.  Alternatively,  numerical  simulations 
of  the  turbulence  structure  provide  PDFs  representative  of  current  conditions 
that  provide  a  means  to  circumvent  the  restrictiveness  of  the  Gaussian 
assumption  and  the  diagnostic  character  of  current  techniques. 

A  third  alternative  is  to  cast  the  Lagrangian  form  (see  equation  (10))  into  a 
differential  equation  that  can  be  solved  as  a  boundary-value  problem. 
Assuming  diffusion  to  be  a  Markov  process  (i.e.,  that  the  current  physical  state 
of  any  particle  is  dependent  only  on  its  immediately  prior  state),  [13]  then 


d<c> 

dt 

+  —(ufat)  <  c(x,t)> )  -  [Kij<  c(x,t)> ] 

ox.  2  oXjdXj 

(11) 

+  S(x,t )  -  r(t)  <  c(x,t)> 

where 

S  = 

a  source  term 

r  = 

a  sink  terms 

Uj  = 

wind  velocity  components 

c  = 

three-dimensional  concentration 

= 

represents  a  diffusion  coefficient. 

We  note  immediately  that  equation  (11)  resembles  the  advective  form  of 
equation  (9)  with  the  diffusion  coefficient  replaced  by  a  tensor  and  subject  to 
the  double  spatial  derivative.  Rigorously,  the  diffusion  coefficient  is  a  tensor. 
The  true  difference  between  Ky  and  D  is  that  Ky  is  dependent  on  Lagrangian 
statistics;  whereas,  D  relies  on  Eulerian  statistics.  The  ability  to  cast  the 
Lagrangian  equation  into  a  boundary-value  problem  does  not  solve  the 
time-space  problem.  The  properties  on  which  Ky  are  dependent  are  difficult  to 
measure  and  are  usually  treated  as  parameters,  the  value  of  which  are  adjusted 
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to  fit  experimental  concentration  measurements.  (The  same  approach  that  is 
used  in  forcing  the  results  from  solutions  to  equation  (10)  to  match 
experimental  data.)  Therefore,  neither  the  equation  nor  the  problem  has  been 
adequately  solved,  and  the  accuracy  of  the  results  are  dependent  on  how  well 
the  physical  conditions  of  concern  correspond  to  conditions  to  which  the  data 
were  fitted. 

An  additional  problem  concerns  solutions  to  equation  (11).  The  diffusion 
equation  is  a  parabolic  differential  equation.  Equations  of  this  type  are 
characterized  by  a  physically  unrealistic  property  that  implies  that  diffusion 
proceeds  at  infinite  velocities,  leading  to  problematic  results  that  must  be 
interpreted.  An  extension  of  the  Markov  process  can  be  used  to  yield  a 
hyperbolic  (telegrapher)  equation  that  describes  diffusion  occurring  at  finite 
velocities.  [14]  For  all  practical  purposes,  the  results  obtained  by  the 
telegrapher  equation  are  minimally  different.  Although  it  solves  one  concern, 
the  hyperbolic  solution  does  not  remove  the  problem  of  defining  K^.  It  was 
suggested  earlier  that  numerical  simulations  provide  a  means  for  obtaining 
applicable  probability  density  functions  for  the  turbulence  (discussion  following 
equation  (9)).  Deardorff  and  Peskin  [15]  and  Deardorff  [16]  demonstrated  a 
modeling  technique  now  called  large-eddy  simulation  (LES)  that  has  the 
capability  to  predict  accurate  turbulence  statistics.  The  model  is  a  solution 
technique  applied  to  the  full  set  of  the  Navier-Stokes  equations,  permitting  the 
full  inclusion  of  meteorology,  terrain,  and  other  factors. 

As  pointed  out  by  Stull,  [17]  the  predictability  achieved  by  any  model  of 
turbulence  structure  must  be  viewed  in  two  ways:  (1)  actual  pattern 
predictability  and  (2)  statistical  functions  from  pattern  forecasts.  While  the 
simulated  patterns  of  the  LES  do  not  necessarily  correspond  to  true  structures, 
it  has  been  shown  that  for  well  designed  simulations,  the  accuracy  of  the 
turbulence  statistics  can  be  maintained  for  forecasts  extending  several  days. 
With  the  ability  to  numerically  produce  a  PDF  for  a  given 
terrain/meteorological  situation,  a  specific  solution  to  either  form  (equations 
(10)  and  (11))  of  the  Lagrangian  formulation  is  possible  without  the  attendant 
ambiguities  introduced  by  parameterization.  The  superiority  of  the  direct 
solution  to  the  diffusion  equation  over  the  Gaussian  plume/puff  technique  was 
demonstrated  by  Lamb  et  al.  [18]  This  report  also  provides  a  short  but 
excellent  theoretical  review  and  discusses  the  unrealistic  Gaussian  assumption. 


6.  Modeling  Transport  and  Diffusion 


There  is  a  very  large  stock  of  hazard  prediction  models  available.  To 
differentiate  among  them,  we  can  separate  the  models  into  three  families: 
(1)  solution  technique,  (2)  level-of-detail  within  the  solution  technique,  and 
(3)  inventory  of  materials  treated.  Solution  techniques  can  be  broadly 
categorized  as:  (la)  statistical,  including  Monte  Carlo  and  particle-in-cell, 
(lb)  parametric,  including  most  Gaussian  puff/plume  models,  or  (lc)  numerical. 
Level-of-detail  considerations  include  such  factors  as  treatment  of  multiple 
sources,  phase  change,  directionality,  meteorology,  and  terrain*.  The  types 
(gas-liquid-solid)  of  material  handled  by  a  model  impact  not  only  the  hazard 
prediction,  because  of  the  range  of  toxicities  and  effects  available,  but  also 
T&D  (particle  mass/density  and  phase  changes  affect  concentrations).  Further, 
because  the  parametric  family  of  T&D  solutions  has  dominated  modeling  for 
military  applications,  discussions  will  focus  on  this  approach.  However,  it  will 
be  immediately  clear  that  the  main  conclusion  is  equally  applicable  to  all 
classes. 


All  hazard  prediction  models  in  use  are  attempts  at  obtaining  solutions  to  the 
advective  form  of  the  basic  diffusion  equation: 
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The  salient  difference  in  solution  techniques  is  the  manner  in  which  the 
stochastic  nature  of  the  atmospheric  processes,  specifically  the  eddy  diffusivity, 
Ky,  is  treated.  It  was  first  recognized  by  Lord  Rayleigh  that  the  Monte  Carlo 
(alternatively  known  as  the  method  of  statistical  trials,  or  the  random  walk 
method)  was  a  suitable  method  for  solving  boundary-value  differential  equations 
of  random  processes.  It  was  not,  however,  until  the  advent  of  high  speed 
computers  that  such  methods  were  realizable.  All  models  in  the  statistical 
family  attack  the  problem  by  examining  the  behavior  of  a  representative  number 
of  entities,  eventually  assigning  the  resulting  average  tendencies  to  the  larger 


*See  appendix  B  for  a  discussion  of  the  interconnection  of  terrain  and  meteorology  as  it 
pertains  to  dispersion  modeling. 
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ensemble  of  particles.  The  distinction  between  the  models  is  in  the  achievement 
of  the  multiple  trial  concept.  Pure  Monte  Carlo  codes  create  the  statistical 
ensemble  through  multiple  runs  on  a  few  particles.  Particle-in-cell  methods 
employ  a  larger  number  of  particles  that  are  initially  assigned  to  cells; 
averaging  is  accomplished  on  the  collection  of  cells.  In  each  case,  the  velocity 
(direction  and  speed)  of  the  individual  particles  during  each  time  step  is 
determined  by  random  functions  (random  number  generators).  To  be  valid,  the 
applications  of  the  random  number  generator  must  result  in  movements  from 
one  location  to  another  that  represent  physical  reality.  To  accomplish  this,  the 
random  number  generator  must  mimic  some  PDF.  The  shape  of  the  probability 
distribution  must  be  empirically  determined  so  as  to  have  the  appropriate 
atmospheric  turbulence  characteristics.  [19]  In  contrast  to  the  random  walk 
techniques,  many  parametric  methods  utilize  particular  solutions  to 
equation  (10).  To  arrive  at  a  particular  solution  a  number  of  simplifying 
assumption  are  necessary. 

A  number  of  simplifying  assumptions  are  necessary  to  determine  a  transition 
probability  for  the  case  of  atmospheric  turbulent  diffusion.  First,  it  must  be 
assumed  that  the  turbulence  is  homogeneous  and  stationary.  In  this  case, 
Q(x,t]x',t')  of  a  particle  depends  only  on  the  displacements  in  time  and  space, 
not  on  when  or  where  the  particle  was  introduced  into  the  flow.  Secondly,  the 
incorrect  assumption  that  the  turbulent  energy  is  normally  distributed  is  applied. 
Following  Monin  and  Yaglom  [12]  a  normal/Gaussian  distribution  of  turbulent 
energy  results  in 

Q(x  -  x',t  -  t)  =  - * — -e<Tp  ‘c/2  (13) 
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where 

=  the  transpose  of  a  column  vector  of  space  correlations  with 
the  components  Q  =  X;-x,'  -  <xj-x/>  (i  =  1,2,3) 

P'1  and  |PJ  =  the  inverse  and  determinant  of  the  matrix  of  which  the 
elements  are  Py  =  <££>. 


An  additional  assumption  of  isotropic  turbulence  results  in  Py  =  0  for  i  *  j.  For 
x  =  t  - 1' 
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Replacing  Q  in  equation  (10).  with  equation  (14),  defining  the  source 
distribution  (a  point  source  is  given  by  the  probability  distribution, 
i) /  =  S5(x-x0)5(t')  dx'),  and  assuming  the  mean  wind  velocity  is  along  a 
coordinate  axis,  we  gain  a  trivariant  joint  normal  time-dependent  solution  of  the 
form  typically  found  in  models  of  the  parametric-analytic  class.  A  specific 
form  (sample  solutions  are  shown  in  appendix  A)  is  dependent  on  the  imposed 
boundary  conditions  as  well  as  other  assumptions.  The  results  for  any  time  t 
are  obtained  through  integration  of  the  appropriate  particular  solution  with 
respect  to  time,  producing  well-behaved  Gaussian  plumes  that  are  seldom 
realized  in  actual  conditions.  The  solutions  are  characterized  by  the  standard 
deviation  of  the  PDF  a  known  also  as  the  Gaussian  spread  factor.  Inspection 
of  the  simple  solutions  (appendix  A)  immediately  reveals  the  causes  for  the 
nonrealistic  results: 

a.  Turbulent  diffusion  (the  entire  point  of  the  model)  is  confined  to  three 
parameters  or  spread  factors  Gr  Additionally,  a  is  a  function  of  time.* 

b.  The  transport  term  U  is  a  constant. 

c.  The  source  S  is  singular  and  constant. 

d.  No  provision  for  loss  of  material  is  included. 

e.  Concentrations  are  relative  to  a  fixed  axis,  rather  than  relative  to  the  center 
of  mass  of  the  puff. 


*The  time  dependence  of  a,  in  reality  a  spatial  and  temporal  dependence,  is  seldom 
treated. 
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Many  different  adjustments  to  the  basic  solutions  have  been  applied,  and  while 
the  purist  may  argue  that  the  solution  is  no  longer  valid;  pragmatically,  such 
results  can  be  forced  to  resemble  experimental  measurements.  To  repair  the 
apparent  violation  of  the  Lagrangian  spatial  correlation  implied  by  e, 
investigators  have  modified  the  equations  to  include  a  Lagrangian  time  scale 
factor.  A  large  variety  of  methods  have  been  applied  to  include  the  effect  of 
material  sinks  (d),  primarily  through  evaporation  and  deposition  as  well  as 
casting  S  as  S(x,t);  thereby,  permitting  time-dependent  additional  sources  (c). 
Such  adjustments  are  valid  because  the  independence  of  x  was  a  requirement 
for  and  time  was  not  a  factor  in  the  derivation  of  equation  (14).  Likewise, 
allowing  U  =  uf(t)  addresses  (b)  but  does  not  invalidate  the  results.  The  most 
severe  deficiency  is  contained  in  a  =  f(Kjj,x,t)  and  the  underlying  assumption 
of  a  normally  distributed  turbulence  field.  Recall,  to  achieve  this  solution,  it 
was  necessary  to  assume  that  the  turbulence  was  h-i-s,  which  for  the  boundary 
layer  is  rarely  a  valid  assumption.  Ignoring  this  problem,  one  is  still  faced  with 
achieving  a  valid  representation  of  turbulence  commonly  through  the  standard 
deviation,  ai5  of  the  normal  (Gaussian)  distribution*.  The  method  most  used 
is  attributed  to  Pasquill,  [4]  and  its  further  development  by  Gifford  and 
others.  [20,21]  These  methods  that  appeal  to  the  physics  embodied  in  the 
Richardson’s  number  and  the  friction  velocity  with  empirical  modifications  for 
the  effects  of  cloud  cover  are  referred  to  as  the  Pasquill-Gifford  (P-G)  stability 
category  methods.  (When  the  downwind  dependence  of  (j;  is  included,  the 
method  is  sometimes  attributed  to  McElroy-Poole.)  [22]  As  the  title  implies, 
these  methodologies  provide  a  means  for  selecting  approximate  spreading 
factors  oi  based  on  estimates  of  the  boundary  layer  stability  condition.  The 
inability  of  the  technique  to  produce  reliable  results  has  led  to  additional,  and 
often  inappropriate,  parameterizations.  Given  a  sufficient  number  of  adjustable 
parameters,  a  model  certainly  should  be  capable  of  reproducing  a  single  given 
result,  providing  grounds  for  claiming  success.  The  prototype  model  suggested 
by  Taylor  and  Pasquill  [4,5,11]  was  designed  for  use  in  the  boundary  layer  as 


‘Atmospheric  turbulence  is  predominantly  non-Gaussian.  For  dispersion,  the  nonuniform 
behavior  is  enhanced  by  the  superposition  of  large  and  small  eddys.  Small  eddys  do  not 
significantly  contribute  to  the  diffusion  process  and  do  not  move  independently  as 
suggested  by  the  Gaussian  assumption.  Rather,  small  eddys  move  within  the  larger 
structures  producing  a  skewed  distribution  of  turbulent  energy. 
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were  parameters  from  mixing-length  theory  that  are  routinely  used  in  attempts 
to  force  a  simple  model  to  replicate  a  complex  process.  Therefore,  if  such  a 
model  is  useful,  its  functionality  must  be  assumed  to  be  only  appropriate  to 
problems  within  the  boundary  layer.  However,  the  boundary  layer  in  general 
and,  specifically,  near  surface  releases  are  the  situations  least  likely  to  encounter 
conditions  adhering  to  the  Gaussian  assumption.  [23]  The  use  of  a  numerical 
methodology  avoids  the  faults  associated  with  the  parametric  approach,  and 
offers  the  possibility  of  a  model  of  greater  general  applicability. 

For  short  periods  of  interest,  some  parametric  approaches  may  provide  an 
adequate  solution  for  uncomplicated  cases  if  the  turbulence  can  be  properly 
characterized.  However,  the  atmosphere  is  a  complex  system  that  eludes  simple 
parameterization,  particularly  for  long  time  periods.  It  was  previously  noted 
that  the  prime  deficiency,  in  the  statistical  and  parametric  approaches,  is  the 
proper  characterization  of  the  turbulence.  With  valid  representations  of 
turbulence  as  a  function  of  time  and  space,  either  equation  (10)  or  (11)  can 
adequately  describe  concentration  of  material  as  a  function  of  time.  It  was  also 
noted  that  properly  executed  LES  models  have  the  proven  capability  of  reliably 
providing  statistics  for  atmospheric  turbulence.  The  LES  is  the  foundation*  of 
a  numerical  approach  to  solving  T&D  problems  and,  ultimately,  the  prediction 
of  hazard  areas  caused  by  the  release  of  toxic  materials  at  any  location,  height, 
or  time.  In  addition  to  correct  statistics  and  removal  of  the  ambiguities  that  are 
inherent  in  parameterizations,  advantages  gained  through  use  of  a  LES  include 
the  following: 

a.  complete  four  dimensional  mean  wind  and  temperature  fields  for  a  more 
complete  treatment  of  transport  and  chemical  reactions 

b.  complete  four  dimensional  treatment  of  terrain  and  diurnal  effects 

c.  prognostic  capability  on  time  and  spatial  scales  pertinent  to  T&D 


*The  LES  is  particularly  suited  for  PBL  applications,  but  proper  structuring  the  code  can 
extend  its  applicability.  This  is  also  true  of  other  methods  for  solving  the  set  of 
meteorological  equations. 
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d.  equal  applicability  to  the  statistical  or  parametric  approach 

In  a  sense,  the  numerical  class  is  not  a  separate  family,  but  a  new  entry  point 
into  the  first  two  families. 


7.  Hazard  Prediction  Modeling,  A  Brief  Survey 


The  point  of  this  report  is  chemical  and  biological  T&D,  the  results  of  which 
are  spatial  concentrations  and  distributions  of  a  given  material.  T&D  is  at  the 
core  of  hazard  predictions,  which  often  leads  to  consideration  of  such  models 
as  atmospheric  dispersion  codes.  Hazard  prediction  generally  requires  a  time 
integration  of  the  concentrations  produced  by  T&D  calculations  yielding 
dosages  that  can  be  assessed  in  terms  of  medical  impacts.  In  this  sense,  T&D 
is  a  module  within  the  hazards  prediction  scheme.  In  assessing  the  applicability 
of  a  given  model  to  a  problem,  this  modular  aspect  should  be  considered. 

The  modular  notion  provides  further  considerations.  Added  modules  will  likely 
increase  computational  time  and  memory  requirements,  but  more  critical  is  the 
proper  and  judicious  combination  of  modules.  To  illustrate  this  point,  consider 
recent  uninformed  attempt  at  increasing  the  reliability  of  hazard  predictions  by 
combining  a  Gaussian  plume/puff  model  with  a  mesoscale  meteorological 
model.  The  Gaussian  model  is  an  analytic  solution  designed  to  provide 
diagnostic  results  for  short  times  and  distances  and  generally  has  a  implicit  grid 
system  fixed  to  the  release  point  with  one  axis  oriented  in  the  direction  of  the 
mean  wind.  The  mesoscale  model  is  a  numerical,  predictive  model,  with  an 
explicit  fixed,  independent  grid.  The  distances  over  which  a  Gaussian  solution 
might  be  considered  valid  are  small  compared  to  the  typical  mesoscale  model 
spacing  of  5  km.  The  entire  domain  of  the  Gaussian  solution  is  likely  to  be 
subgrid  scale  and  unresolvable,  relative  to  the  mesoscale  meteorological  model. 
More  important,  and  increasingly  so  when  higher  resolution  meteorological 
models  are  used,  is  the  incompatibility  of  the  methodology  used  to  describe 
turbulence,  the  primary  influence  on  changing  concentrations.  The 
incompatibility  results  in  neglect  of  the  available  turbulence  information  or  the 
need  for  substantial  additional  computations,  which  defeats  the  simplicity  of  the 
analytic  model.  Failure  to  use  the  turbulence  information  negates  any 
advantage  anticipated  from  the  additional  overhead  associated  with  the 
improved  meteorology.  The  complexity  involved  in  using  the  data  suggests 
replacing  the  Gaussian  technique,  which  is  not  always  applicable,  with  a 
compatible  numerical  solution  of  the  basic  diffusion  equation.  The  value  of 
improved  meteorology  can  then  be  scenario  driven. 
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A  point  of  particular  relevance  for  military  application  is  the  persistency  of 
agents.  As  toxic  agents  become  more  persistent,  and  delivery/intercept 
mechanisms  become  ubiquitous,  the  need  for  models  capable  of  adequate 
treatment  of  transport  and  diffusion  for  long  time  periods  and  at  all  levels  in  the 
atmosphere  become  more  important.  For  releases  of  persistent  agents  at 
mid-to-high  altitude,  prognostic  or  predictive  models  are  increasingly  desirable. 
The  following  review  includes  models  that  have  been  used  by  various  military 
agencies,  but  is  not  an  exhaustive  compilation.  However,  after  the  family  to 
which  a  model  belongs  can  be  identified,  a  cognizant  user  can  immediately 
determine  the  level  of  accuracy  to  be  expected  for  a  specific  application.  For 
additional  details  on  a  specific  model,  the  reader  is  directed  to  cited  references. 

7.1  Statistical  Models 

Models  in  the  statistical  (Monte  Carlo)  family  (MESO,  [24]  RAPTAD,  [25] 
and,  ADPIC  [26])  have  general  applicability.  As  pointed  out,  the  accuracy  of 
statistical  T&D  results  is  dependent  primarily  on  the  accurate  portrayal  of  the 
turbulence  data.  Capable  of  diagnostic  results,  these  models  are  well  suited  to 
running  in  a  predictive  mode  if  coupled  to  a  prognostic  meteorological  model. 
RAPTAD  is  the  only  statistically  based  T&D  model  linked  with  such  a  model. 
However,  in  general  terms  the  associate  model  HOTMAC  is  capable  of  a 
minimum  grid  resolution  of  2  km*,  which  is  questionably  too  large  for 
appropriate  turbulence  statistics.  Additionally,  RAPTAD  does  not  address 
issues  associated  with  change  of  phase;  therefore,  it  is  restricted  to  use  for  the 
dispersion  of  gaseous  material.  The  Monte  Carlo  approach  is  a  viable  and 
attractive  method  for  the  solution  of  the  hazard  prediction  problem.  Its  primary 
drawback  is  the  fact  that  its  accuracy  goes  as  NJ/\  This  suggests  that  a  very 
large  number  of  statistical  trials  N  must  be  computed.  Each  of  these  models  are 
dispersion  rather  than  hazard  prediction  codes. 


‘General  applicability  demands  that  a  model  treat  the  majority  of  expected  conditions. 
Meteorological  models  developed  with  a  "hydrostatic"  assumption  (e.g.,  HOTMAC)  suffer 
from  the  inability  to  treat  convectively  unstable  atmospheric  conditions.  For  such 
conditions,  it  is  not  advisable  to  use  a  computational  grid  with  spacing  smaller  than 
approximately  2  km. 
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7.2  Gaussian  Models 


The  comparison  of  models  within  the  parametric  family,  including  OB/DG 
(Ocean  Breeze/Dry  Gulch),  CHARM  (Chemical  Agent  Release  Model),  AFTOX 
(Air  Force  Toxic  Corridor),  NUSSE  (Non-Uniform  Simple  Surface 
Evaporation),  GAPCAP  (Generation  of  Assessment  Patterns  for  Clouds  of 
Airborne  Particles),  VLSTRACK  (Vapor  Liquid  Solid  Tracking),  MACH-LT 
(Model  of  Atmospheric  Chemical  Hazards  -  Laptop) ,  requires  a  comparison  of 
the  level  of  detail  included  in  each  of  the  models.  (See  appendix  B  for  a 
discussion  of  terrain  interaction  in  T&D  modeling.)  All  the  mentioned  models, 
except  OB/DG,  rely  on  the  P-G  stability  parameterization  to  determine  the 
effective  atmospheric  diffusion.  The  utilization  P-G  categories  limits  the 
applicability  of  a  model  to  descriptions  of  releases  on  or  near  the  surface  of  the 
earth. 

1.  OB/DG  -  [27]  Look-up  tables  determine  atmospheric  effects;  tables  available 
for  36  gaseous  materials;  gases  considered  are  neutrally  buoyant.  Continuous 
point  sources  only.  No  mass  loss  or  gain  considered.  Single  point  and  time 
wind  and  temperature  input.  Flat  terrain.  Diffusion  parameterized  to  single 
vertical  temperature  variation  with  height.  Dispersion  code. 

2.  CHARM  -  [28]  Standard  Gaussian  Puff  routine  with  parameterized  increase 
in  dispersion  in  wake  of  obstructions  (buildings).  Gases  only  -  data  available 
for  136  gases  with  wide  variety  of  gas  densities.  Continuous  point  sources 
only.  Mass  gain  from  evaporation  allowed.  Single-point  wind  and  temperature 
at  multiple  times.  Flat  terrain.  Diffusion  based  on  criteria  similar  to  P-G 
methodology.  Dispersion/hazard  code. 

3.  AFTOX  -  [29]  Standard  Gaussian  Puff  routine.  Handles  all  chemicals  with 
known  molecular  weight.  Gases  only,  but  treats  mass  gain  from  evaporation. 
Single  point  and  time  wind  and  temperature.  Flat  terrain.  Diffusion  based  on 
P-G  methodology,  boundary-layer  parameterization  used  to  calculate  stability. 
A  dispersion  code. 

4.  NUSSE4  -  [30]  Standard  Gaussian  Puff  routine.  Continuous-point  sources 
only.  Tracks  liquid  and  gaseous  clouds.  Losses  caused  by  deposition  and 
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(vapor)  mass  gain  from  surface  evaporation  of  deposited  liquids  treated. 
Extensive  use  of  agent  look-up  tables.  Single-point  wind  and  temperature.  Flat 
terrain.  Diffusion  based  on  P-G  parameterization.  Dispersion/hazards  code. 

5.  GAPCAP  -  [31]  This  model  is  a  variation  of  the  NUSSE  family  of  models, 
modified  to  treat  solid  biological  agents,  accounting  for  the  associated  decay  of 
such  agents.  Uses  look-up  tables  for  agent  parameters.  Flat  terrain. 
Single-point  wind  and  temperature.  Diffusion  based  on  P-G  parameterization. 
Dispersion/hazards  code. 

6.  VANDECAP/PACCE  -  [32]  Modified  version  of  GAPCAP.  Flat  terrain. 
Attempts  to  use.  a  single  profile  of  wind  and  temperature  to  permit  application 
to  high  altitude  releases.  Dispersion/hazards  code. 

7.  VLS TRACK  -  [33,34]  Standard  Gaussian  Plume/Puff  routine. 

Continuous-point  and  line  sources.  Tracks  clouds  consisting  of  liquid,  gaseous, 
and  solid  phases.  Losses  caused  by  deposition  and  mass  gain  from  evaporation 
surface  evaporation  of  deposited  liquids  treated  in  some  versions.  Agents  in 
database.  Various  input  schemes  for  wind  and  temperature  (discussed  below). 
Flat  terrain.  Dispersion  based  on  P-G  parameterization.  Dispersion/hazards 
code. 

8.  MACH-LT  -  [35]  Modified  Gaussian  Plume/Puff  routine.  Continuous-point 
and  line  sources.  Tracks  liquid  and  gaseous  and  solid  clouds.  Losses  caused 
by  deposition  and  mass  gain  from  evaporation  surface  evaporation  of  deposited 
liquids  treated.  Agents  in  database.  Single-point  wind  and  temperature.  Flat 
terrain.  Diffusion  based  on  statistics  derived  from  LESs  over  similar  terrain. 
Diffusion/hazards  code. 


7.3  Numerical  Models 


Despite  the  astounding  advances  in  computation  power  experienced  over  the 
past  decade,  a  complete  numerical  description  of  all  the  physical  processes  in 
a  complete  first  principles  approach  is  not  yet  possible  for  the  scales  that  would 
be  required  for  very  accurate  or  reasonable  hazard  prediction  modeling. 
However,  it  is  possible  to  use  numerical  techniques  that  provide  a  substantial 
improvement  to  either  the  representation  of  the  problem  or  the  application  of 
the  statistical  techniques. 

1.  SCIPUFF  -  [36]  Second-order  Closure  Integrated  PUFF.  Continuous-point 
and  line  sources.  Tracks  clouds  consisting  of  liquid,  gaseous,  and  solid  phases. 
Treats  mass  losses  and  gains.  A  direct  numerical  solution  of  the  diffusion 
equation  (12)  with  turbulence  treated  in  a  manner  analogous  to  standard 
perturbation  theory.  Capable  of  calculations  for  complex  terrain.  Diffusion 
is  based  on  statistics  for  turbulence  spectra  inferred  from  climatological  data 
and  more  recently  derived  from  a  turbulence  closure  model.  The  primary 
limitation  of  SCIPUFF  is  in  the  treatment  of  K  the  diffusion  coefficient  as  a 
scaler  quantity.  Dispersion  code. 

2.  ABCSIM  -  [35]  Atmospheric  Biological  Chemical  Simulations. 

Continuous-point  and  line  sources.  Tracks  clouds  consisting  of  liquid,  gaseous, 
and  solid  phases.  Treats  mass  losses  and  gains.  A  family  of  modules  that 
include  an  LES  complex  terrain  meteorological  model  and  a  direct  numerical 
solution  of  the  diffusion  equation  (12).  Turbulence  statistics  are  derived  from 
the  LES,  and  the  diffusion  coefficient  is  treated  in  its  tensor  form. 
Diffusion/hazards  code. 

3.  ABCSIM-LT  -  [35]  Atmospheric  Biological  Chemical  Simulations  -  Laptop. 
A  scaled-down  version  of  ABCSIM  for  use  on  small  computers.  Reduced 
capability  relates  to  the  computation  of  the  vertical  structure  of  wind  and 
turbulence.  Speed  and  memory  requirements  restrict  LES  computations  to 
approximately  the  lowest  1  km  of  the  atmosphere.  Transport  and  diffusion 
calculations  are  unchanged.  Diffusion/hazards  code. 
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8.  Conclusions  and  Recommendations 

The  inference  of  this  report  is  that  adequate  predictions  of  hazards  associated 
with  the  dispersion  of  agents  into  the  atmosphere  requires  that  substantial 
attention  be  placed  on  the  efficacy  of  the  method  in  which  T&D  are  treated. 
While  VLSTRACK  is  not  a  T&D  model,  it  does  contain  and  is  under 
consideration  to  be  the  model  of  choice  by  many  DoD  agencies.  Therefore,  a 
number  of  specific  comments  are  addressed  to  that  model. 

A  description  of  the  evolution  of  a  cloud  of  material  suspended  in  the 
atmosphere  can  be  separated  into  transport  (the  movement  of  the  center  of 
mass)  and  diffusion  or  the  changing  distribution  of  the  material  about  the 
centroid.  Transport  is  derived  from  knowledge  of  the  wind  at  the  location  of 
the  centroid  at  given  times,  or  estimated  by  assuming  a  mean  wind  over  the 
entire  travel  distance.  Either  case  is  the  result  of  a  relatively  simple  calculation. 
Superimposed  on  this  are  other  physical  processes  such  as  phase  change  and  the 
more  elusive  process  of  turbulence-induced  dispersion. 

In  general,  diffusion  models  based  on  a  parametric  paradigm  are  useful  as 
estimators  or  indicators  of  trends,  provided  atmospheric  conditions,  including 
terrain  effects,  are  relatively  benign  and  constant  and  dispersion  times  are 
relatively  short.  Persistent  agents  and  complex  terrain  (appendix  B)  can  only 
be  addressed  at  increasingly  lower  accuracy  levels  and  at  considerable 
computational  expense.  The  Gaussian  model,  while  physically  reasonable  for 
pure  molecular  diffusion,  is  not  generally  representative  of  boundary-layer 
turbulent  diffusion  except  for  large  averaging  times.  [23]  Therefore, 
Gaussian-based  models  have  a  severe  inconsistency  between  the  underlying 
physical  realities  and  modeling  assumptions.  The  basis  of  the  model,  the 
diffusion  equation  (12),  is  derived  from  analyses  of  mo.lecular  motion  and,  as 
such,  is  dependent  solely  on  physical  properties  of  the  fluid.  Turbulent 
diffusion  is  not  intrinsic  to  the  fluid  or  the  dispersing  species,  rather,  it  is  driven 
by  external  forces.  Without  regard  to  the  questionable  assumptions  required  to 
obtain  the  Gaussian  solution,  this  suggests  caution  should  be  exercised  in  the 
application  of  the  equation  for  the  description  of  events  dominated  by  the 
mechanism  termed  turbulence.  The  limited  success  of  early  investigators  who 
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attempted  to  extended  the  notions  of  molecular  diffusion  to  include  turbulence 
effects  masks  the  fragility  of  a  generalization  of  the  approach. 

In  addition,  diffusion  resulting  from  turbulence  is  a  local  phenomenon,  and 
unless  the  cause  of  the  turbulence  is  uniform  over  a  large  area  (or  minimally 
the  path  on  which  the  calculation  is  occurring),  utilizing  a  parameterization 
based  on  measurements  at  a  single  location  is  subject  to  large  errors. 
Measurements  used  to  verify  models,  including  early  models  based  on  the  P-G 
stability  category,  were  conducted  in  regimes  that  adhered  to  this  limitation. 
It  is  not  surprising  that  there  is  some  agreement  between  model  results  and 
measurements.  These  comparisons  also  verify  the  decrease  in  accuracy  with 
time.  This  is  an  expected  result,  a  single  simple  solution  can  not  be  expected 
to  adequately  track  the  effects  of  a  dynamic  atmosphere  and  Gaussian 
plume/puff  models  are  based  on  one  simplified  solution  of  the  molecular 
diffusion  equation.  Coupling  such  solution  methods  to  a  four-dimensional 
dynamic  meteorological  model  has  questionable  merit. 

Meteorological  models  useful  for  application  to  T&D  problems  add  substantial 
computational  overhead.  If  there  is  a  benefit  to  be  gained  from  the  additional 
computations,  it  should  be  in  a  better  representation  of  the  mechanism  that 
contributes  to  the  dilution  of  the  species  being  tracked.  If  the  meteorological 
model  is  used  solely  for  wind  information,  the  more  important  information  is 
being  discarded.  Using  the  available  turbulence  information  to  improve  a 
Gaussian  model  would  require  considerable  restructuring,  and  the  benefit  is 
questionable.  Certainly,  a  better  estimate  of  the  standard  deviation  a  should 
improve  model  specificity  and  results,  but  the  model  remains  Gaussian  and  not 
representative.  The  additional  computations  and  necessary  factor  fiddling 
contradicts  the  purpose  of  an  analytic  solution.  The  contradictions  are 
intensified  when  a  meteorological  model  of  sufficient  resolution  is  used. 

Consideration  of  resolution  is  an  important  aspect  because  it  controls 
computational  resource  requirements  as  well  as  accuracy  of  results.  For 
example,  consider  a  case  used  to  validate  VLSTRACK.  The  largest  total  area 
coverage  for  any  trial  was  about  41,000  m2,  equivalent  to  a  1 14-m-radius  circle. 
Using  the  HOTMAC  mesoscale  model  with  an  approximate  2000-m  minimum 
resolution  (refer  to  footnote  on  a  page  28),  a  single  HOTMAC  cell  is  then 


97  times  larger  than  the  largest  area  used  in  the  validation  of  the  VLSTRACK 
model.  Therefore,  at  best,  information  from  two  HOTMAC  grid  points  would 
be  useful  to  the  entire  T&D  calculation.  Because  VLSTRACK  is  not 
configured  to  use  turbulence  information  from  a  meteorological  model,  only 
wind  information,  approximately  20  percent  of  the  relevant  information,  would 
be  used.  From  a  diagnostic  point  of  view,  the  single  data  point  input  method 
of  the  NUSSE  model  would  have  been  sufficient.  HOTMAC  (or  any  mesoscale 
model)  used  in  conjunction  with  VLSTRACK  can  provide  prognostic  winds  but 
at  a  resolution  insufficient  to  substantially  influence  the  results.  Models  for 
higher  resolution  winds  are  available,  but  the  computational  overhead  required 
to  use  these  products  with  VLSTRACK  is  inefficient. 

Additionally,  users  of  the  model  are  cautioned  against  reliance  on  the  validation 
and  documentation  of  VLSTRACK.  [33,34]  These  works  were  incorrectly 
titled  validation.  The  work  reported  in  these  nonreferred  records  documented 
a  parameter-fitting  methodology.  The  VLSTRACK  validation  reported  by 
Sterle  et  al.  [37]  documents  the  inability  of  the  analytic  Gaussian  model, 
specifically  VLSTRACK,  to  replicate  field-trial  data  collected  during  less  than 
stressful  terrain  and  meteorological  conditions.  Sterle  et  al.  attributes  the  less 
than  acceptable  performance  of  the  model  to  "inherently  large  variability  and 
unpredictability  of  the  meteorological  conditions,  which  play  a  critical  role  in 
determining  the  distribution  pattern  of  agents  and  stimulants  following 
dissemination  into  the  atmosphere.  Users  of  VLSTRACK  or  any  other 
chemical/biological  dispersion  model  should  be  cognizant  of  this  limitation  and 
thereby  avoid  placing  too  much  weight  on  the  model  outputs  most  affecting 
borderline  decision  making."  [37] 

It  must  be  pointed  out  that  such  cautions  should  be  applied  in  any  case,  but 
meteorology  is  not  completely  unpredictable  for  the  time  scales  of  interest,  and 
the  cited  limitations  are  applicable  only  to  highly  parametric  Gaussian  models 
in  which  all  the  variability  and  unpredictability  of  the  meteorological  conditions 
is  contained  in  a  single  parameter.  Likewise,  the  abandonment  of  models  such 
as  VLSTRACK  would  focus  funding  and  labor  on  more  reasoned  approaches 
to  the  proper  inclusion  of  variable  meteorology  and  turbulence.  Regardless  of 
the  paradigm  used,  an  adequate  description  of  the  diffusion  process  can  only  be 
obtained  if  the  distribution  of  turbulent  energy  is  known  or  correctly  estimated. 
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Models  based  on  the  Gaussian  assumption  are  confined  to  a  normal  distribution 
and  to  isotropicity,  homogeneity,  and  stationarity,  which  do  not  have  general 
applicability  in  the  boundary  layer  of  the  atmosphere.  In  the  absence  of  proper 
turbulence  statistics,  statistical  and  numerical  models  can  use  a  Gaussian 
assumption  but  are  normally  free  to  apply  any  distribution  to  the  turbulent 
energy.  The  accuracy  of  these  simulations  is  dependent  on  the  number  of  trials 
and/or  the  shape  of  the  energy  distribution.  In  this  aspect,  parametric  models 
are  the  least  flexible  and  subject  to  the  largest  errors. 
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Acronyms  and  Abbreviations 


ABCSIM 

Atmospheric  Biological  Chemical  Simulation 

ABCSIM-LT 

Atmospheric  Biological  Chemical  Simulation  -  Laptop 

AFTOX 

Air  Force  Toxic  Corridor 

CHARM 

Chemical  Agent  Release  Model 

DoD 

Department  of  Defense 

GAPCAP 

Generation  of  Assessment  Patterns  for  Clouds  of  Airborne 
Particles 

h-i-s 

homogeneous-isotropic-stationary 

LES 

large-eddy  simulation 

MACH-LT 

Model  of  Atmospheric  Chemical  Hazards  -  Laptop 

NUSSE 

Non-Uniform  Simple  Surface  Evaporation 

OB/DG 

Ocean  Breeze/Dry  Gulch 

PBL 

planetary  boundary  layer 

PDF 

probability  density  function 

P-G 

Pasquill-Gifford 

Re 

Reynolds  number 

SBL 

surface  boundary  layer 

SCIPUFF 

Second-order  Closure  Integrated  PUFF 

T&D 

transport  and  diffusion 

VLSTRACK 

Vapor  Liquid  Solid  Tracking 
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Appendix  A 


Gaussian  Solutions  for  Lagrangian  Dispersion 


The  following  are  examples  of  solutions  to  the  diffusion  equation  (12)  using 
Gaussian  assumptions  of  equation  (14)  for  mean  concentrations: 

1.  Instantaneous  point  source  at  Xq.  Mean  wind  %  -  U;  source  strength, 
S  (g),  no  boundary  conditions. 

(*  -  %  -  -  V*  C*3  -*3p)2 

«**,*)> - - /  2o'J(t)  20j2w  2<”2(,)  (A-1) 

(27t)3o1(f)o2(f)o3(r) 

2.  Continuous-point  source  at  Xq.  Mean  wind  Uj  =  U;  source  strength,  S  (g/s), 
no  boundary  conditions. 

(*2  -*2/  _  ~  *3/ 

S  W<*  -  V  -  W  (A-2) 

<c(x.t)  >  -  — - £ 

(2tc)o2(x  -  xlo))o3(x,  -  x{)U 


3.  Instantaneous  point  source  at  Xq.  Mean  wind  u;  =  U;  source  strength, 
S  (g),  reflecting  barrier  (ground)  at  x3  =  0. 


(x  -  -  utf 


<c(x,t )  >  = 


2oj  (f) 


(2Tc)3a1(0a2(0a3(0 


(*2  -  X^)2  (*,  -  X3q)2  _  +  V2 

2a22(0  ][e~  W(t)  +  ^  2^  -V] 


(A-3) 


4.  Continuous-point  source  at  Xq.  Mean  wind  uf  —  U;  source  strength,  S  (g/s), 
reflecting  barrier  (ground)  at  x3  =  0. 


_ s _ 

(2 n)a2(x  -  x,o))o3(Xj  -  x^U 


(*2  -  V2  __  <*3  ~  *3/ 

2oI2C*  -  *lo)j.g  2<T32(3t,  -  X1(>) 


(*3  - 

2o32(xi  -  x,0)_j  (A-4) 


<c(x,r)  >  = 


c 


+  e 


Appendix  B 

Aspects  of  Terrain  Influence  on  Dispersion  Modeling 
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The  dispersion  of  material  release  near  the  ground  is  subject  to  a  number  of 
factors  of  which  all  but  one  are  correlated  to  meteorological  considerations. 
The  exceptional  influence  is  the  fact  that  the  surface  provides  a  no  flux 
boundary  condition.  In  numerical  and  statistical  models,  this  condition  is  easily 
treated  in  the  natural  grid  system.  In  the  case  of  the  parametric  models, 
equations  (A-3)  and  (A-4)  of  appendix  A  illustrate  the  methodology.  Except 
in  flat  terrain,  this  procedure  adds  an  additional  complexity  that  is  not  easily 
treated. 

Terrain  influences  on  meteorology  that  are  of  major  significance  to  dispersion 
are  (a)  steering  effect,  (b)  generation  of  mechanical  turbulence,  (c)  thermal 
effects  on  the  general  flow,  the  turbulence,  and  possibly  the  chemistry.  The 
steering  influence  of  the  terrain  on  the  mean  air  flow  is  the  easiest  to  imagine; 
and,  as  the  discussion  in  the  main  text  suggests,  the  simplest  to  include  in  any 
type  of  dispersion  model.  The  parametric  class  simplicity  demands  an 
assumption  of  constant  wind  direction.  However,  the  nonconstancy  of  wind 
direction  is  a  major  reason  for  the  inclusion  of  terrain  effects.  Without  this 
assumption,  the  same  invalidation  or  stricture  associated  with  the  inclusion  of 
surface  roughness  effects,  discussed  next,  are  applicable. 

The  generation  of  mechanical  turbulence  is  a  function  of  the  surface  roughness. 
Thermally  generated  turbulence  is  dependent  on  surface  composition.  Both  are, 
in  general,  localized  effects  and  a  most  cogent  reason  for  including  terrain 
effects  in  a  dispersion  model.  The  simplicity  of  a  parametric  model  is 
destroyed  when  attempts  are  made  to  include  the  general  variability  associated 
with  terrain  roughness  and  directional  variability.  In  addition  to  the 
requirement  to  use  an  alien  grid  system  which  substantially  increases  the 
necessary  number  of  computations,  either  the  validity  of  the  equations  (the  form 
in  general  use  assumes  that  the  wind  direction  is  constant  and  along  one  of  the 
solution  axis)  is  destroyed  or  the  need  to  use  a  more  general  and  complicated 
solution  to  equation  (12)  is  dictated.  Current  parametric  models  are  not 
designed  to  incorporate  the  variability  of  turbulence  effects.  All  models  of  this 
class  negate  the  advantages  of  including  terrain  through  reliance  on  the 
generalization  represented  in  the  Pasquill-Gifford  stability  categories. 


Thermal  effects  associated  with  diurnal  features  such  as  mountain/valley  winds 
and  land/sea  breezes  are  additionally  compelling  and  important  reasons  for  the 
inclusion  of  terrain  effects. 

The  foregoing  demonstrates  that  the  terrain  only  indirectly,  through  the 
meteorology,  influences  transport  and  diffusion.  Therefore,  with  the  single 
exception  of  providing  a  reflecting  barrier  in  the  calculation  of  vertical  transport 
and  diffusion,  terrain  is  the  providence  of  the  meteorological  model. 
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Appendix  C 


Comments  on  the  Model  VLSTRACK 


The  Vapor  Liquid  Solid  Tracking  (VLSTRACK)  documentation  has  numerous 
errors  and  inaccuracies,  that  would  and  could  not  escape  the  attention  of  a 
serious  scientific  reviewer.  The  document  is  useful  as  a  user’s  manual,  but  the 
scientific  content  reveals  a  lack  of  complete  understanding  of  meteorology  and 
of  T&D  theory.  The  procedures  for  model  validation  as  exhibited  in  the 
validation  document  does  not  meet  the  requirements  of  AR-5-11,  and  further 
demonstrates  a  lack  of  familiarity  with  the  difference  between  model  validation 
and  model  testing.  It  is  unlikely  that  the  model  is  capable  of  results  sufficiently 
better  than  non-uniform  simple  surface  evaporation  (NUSSE4).  The  published 
differences  in  results  [33,34]  are  most  likely  due  to  the  manipulation  of 
VLSTRACK  modeling  factors,  during  "validation",  leading  to  a  minimization 
of  the  difference  between  field  data  and  model  results.  No  such  manipulations 
were  performed  during  the  NUSSE4  model  evaluations  .  In  this  sense  as  well 
as  in  a  more  specific  sense,  versions  of  VLSTRACK  later  than  1.2  have  not 
been  validated.  In  particular,  validation  of  versions  that  contain  significant 
changes  to  accommodate  the  higher  resolution  meteorology  has  not  been 
accomplished.  The  types  of  improvements  necessary  to  gain  an  advantage  from 
"better/representative"  meteorology  can  be  more  efficiently  applied  to  either  the 
statistical  or  numerical  models.  Further,  if  as  noted  in  Sterle  et  al.,  [37]  no 
significant  differences  exist  in  the  results  of  version  1.4  and  2.0,  one  must 
question  the  need  for  the  "improvements"  as  well  as  the  significance  of 
coupling  VLSTRACK  and  a  meteorological  model,  and  the  increased 
complexity  such  a  combination  introduces  with  no  change  in  results. 

It  is  recommended  that  if  VLSTRACK  is  used,  that  it  is  used  in  a  stand  alone 
configuration  (VLSTRACK  1.2),  and  that  no  further  development  be  applied. 
The  model  is  state-of-the-art  for  its  configuration  and  "improvements"  can  only 
provide  minimal  changes. 

In  assessing  the  value  of  a  model,  one  must  be  acutely  aware  of  its  limitations 
and  strengths.  Hazard  prediction  simulations  are,  by  nature,  multi-disciplined, 
requiring  elements  of  physical  chemistry,  meteorology,  atmospheric  physics, 


Tn  comparison  to  NUSSE,  VLSTRACK  does  provide  a  easier-to-use  user  interface  for 
model  initialization  and  it’s  graphical  output  substantially  eases  the  task  of  output 
analyses. 
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medicine  and  some  esoteric  science  that  is  capable  of  describing  the  initial 
conditions.  To  expect  a  single  individual  to  possess  the  all  the  required 
knowledge  to  produce  an  adequate  model  for  the  protection  of  personnel  is 
courting  disaster. 
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